High intensity focused ultrasound (HIFU) is widely applied as a non-invasive technique in clinical diagnose and therapy, including hemostasis, thrombolysis, shock wave (SW) treatments, ultrasonic drug delivery, and lately in treatment for painful conditions.^[1–5] Taking advantage of the non-invasive, small, well-circumscribed thermal lesion, extracorporeal shock wave lithotripsy (ESWL) has been widely applied for stone disease treatments in the past few decades.^[3,6–8] Despite these impressive therapeutic effects, there is non-ignorable evidence that ESWL causes tissue injury, such as hematomas, hematuria, renal and perirenal hemorrhage, and kidney enlargement.^[9,10] Including the above tissue damage caused by cavitation, bubbles formed in the propagation path scatter and absorb much of the acoustical energy before it reaches the target area. Furthermore, it also changes the focal zone location in the target tissue.^[11] Many efforts have been made to control the cavitation in SW-induced treatments.^[6,12–15] It is reported that the bubble cavitation could be suppressed by a small overpressure or slow pulse repetition frequency during ESWL.^[12,13] Evan et al. pointed out that kidney damage and renal functional changes could be minimized by waveform control that reduces cavitation in shock wave lithotripsy (SWL).^[6] Zhong and Zhou developed an in-situ pulse superposition technique to suppress large intraluminal bubble expansion without weakening the curative effect.^[14,15]

Although primary cavitation excited by SWs is crucial to the stone comminution, the residual daughter bubbles accompanying the collapse of cavitation bubbles also play an important role.^[7,16–19] These residual daughter bubbles can persist from one shock to the next, in the order of one second.^[18] On one hand, these bubbles will attenuate the amplitude of the next SW and delay the time for the SW to reach the stone,^[13,19] which leads to weakening stone comminution efficiency. The influence is more significant at higher pulse-repetition frequency (PRF), as less time is available for residual daughter bubbles to passively dissolve between sequent SWs. In order to enhance the comminution efficiency, Duryea et al. developed an active residual bubble removal method.^[20,21] They suggested that adjusting the cavitation environment can improve comminution efficiency when PRF is at a high rate. On the other hand, residual bubbles may also act as nuclei that can be excited by subsequent pulses. Therefore, the residual bubbles may improve the therapeutic effect,^[7] or cause further tissue damage. Therefore, it may be favorable to achieve an optimal SW treatment^{[7,13, 22]} by adjusting the dissolving process of residual cavitation bubbles between two sequent SWs.

It is known that free gas bubbles present in a liquid are usually unstable, which means they will tend to dissolve due to the external pressure generated by surface tension that will force the gas out of the bubble into solution in the liquid. Epstein and Plesset developed equations that can be used to calculate the expected lifetime of an air bubble in an unbounded liquid.^[23] Effects have been done to observe the UCA dissolution process by the ACD system during in vitro experiments.^[24] However, the dissolution process of gas bubbles generated in vivo is not yet well known. Therefore, we intended to investigate the dissolution processes of cavitation bubbles generated during in vivo SW treatments. The general hypotheses tested in this study were: (i) B-mode imaging is an effective technique to monitor the dissolution procedure of cavitation bubbles; and (ii) the cavitation bubble dissolution procedure detected by B-mode ultrasound can be simulated based on the calculation of scattering cross-section with employing gas bubble dissolution equations and Gaussian bubble size distribution. Firstly, in vitro experiments were conducted to monitor the dissolution process of biSpheres contrast agent bubbles using both the single-transducer ACD system and B-mode ultrasound imaging system. To verify the efficiency of the B-mode imaging system, its collected signals were compared with the signals detected by the ACD system, which had been proved to be capable of detecting microbubbles destruction and the subsequent dissolution process.^[24] Furthermore, the normalized decay curves of hyperechoic regions measured by the B-mode imaging system^[25] during both in vivo ESWL treatments of pig kidneys and in vivo ESWT treatments of plantar fasciitis. The results showed that monitoring bubble dissolution procedures using B-mode imaging could provide a better understanding of the temporal evolution of cavitation bubbles after SW pulses ceased in vivo. In addition, the characteristics of cavitation bubbles, such as the bubble size distribution and gas diffusion, can be estimated by simulating the experimental data properly.

The purpose of this paper was to simulate the cavitation bubble dissolution process detected after SW treatments in ESWL and ESWT. Because the experimental data were collected using B-mode ultrasound, the simulation should be based on the calculation of the backscattering cross-section. Furthermore, the simulation should also include the effects of bubble size distribution, because the detected signals involved the scattering from a distribution of bubbles. Therefore, the final simulation should incorporate the changing radius of the dissolving bubbles, R(t), into the scattering cross-section equations, taking into account the bubble size distribution. Chen et al. ^[24] performed just this type of analysis using standard ACD in an in vitro setting by monitoring the destruction of ultrasound contrast agents (UCA) microbubbles. The results were quite good, and thus we will use their simulation and analysis techniques for our purposes.

2.1. Bubble dissolution

It is known in an in vivo setting, bubble dissolution is complicated because they can be stabilized by other surfaces or coatings, or they might move because of blood flow, or there might be bubble-bubble interactions. We will simplify the process by assuming that the dissolution is from an unbounded liquid, without interactions. We hypothesize that discrepancies between the simulations and data are most probably due to the effects mentioned. The dissolution rate of a gas bubble in liquid can be calculated by the following equation:^[23]

(1)

where τ = 2Mσ/R _G T, t is time, C _∞ is the concentration of dissolved gas in solution, C ₀ is the saturation concentration of dissolved gas, and _M is the molecular weight of gas.

Other parameters are defined in Table 1. The initial bubble radius condition used in the calculation was R ₀, which corresponded to the equilibrium bubble radius at t = 0.

Table 1.

Table 1.

Table 1. Notations used for bubble dissolution calculation. . Variable Parameter Simulation value D A diffusion constant of air 2×10–9 m2/s L A Ostwald coefficient of air 1.71×10–2 R G gas constant 8.314 J/mole·K T temperature 300 K σ surface tension 0.07275 N/m ρ g gas density 1.161 kg/m3 η L liquid shear viscosity, from Church et al., 1995[26] 0.001 Pa·s ρ L density of a Newtonian liquid 103 kg/m C p heat capacity at constant pressure, from CRC handbook 240.67 K g thermal conductivity, from CRC handbook, for air at 300 K and 1 atm 0.00626 c acoustic velocity in the liquid 1500 m/s γ gas adiabatic constant 1.07

Table 1. Notations used for bubble dissolution calculation. .

2.2. Calculation of bubble scattering cross-section

Once we know the size of a bubble at any time, we can then calculate its scattering cross section as a function of time based on the assumption of the scattering from small particles with the wavelength λ ≫ R:^[27]

Thermal damping constant:

Radiation resistance damping constant

Viscosity damping constant

(7)

where a _e is the equilibrium bubble radius, ω is the frequency, ω ₀ is the resonance frequency, and P ₀ is the hydrostatic pressure. Other parameters are listed in Table 1.

With a scaling factor, the square root of the scattering cross section can be used to simulate the expected scattering amplitude vs. time that we would see on a B-mode image.

2.3. Size distribution of bubbles

During ESWL or ESWT treatment, bubble ‘clusters’ are generated by shock wave pulses, rather than single bubbles. The size distribution of these bubbles is not homogeneous. For simplicity, a Gaussian function^[24] with mean radius of R ₀ and standard deviation of SD will be employed to cavitation bubble size distribution W (R ₀)

(8)

where RR is in range of radii.

The bubble size distribution was considered as a weighting factor when the scattering cross section was calculated. Thus, the final simulated scattering cross-section for a distribution of bubbles is given by:

It was proved that the single-transducer ACD system was capable of detecting UCA microbubble destruction and the subsequent dissolution process.^[24] Here, the dissolution of biSpheres bubbles was observed in vitro using both ACD and B-mode imaging. The comparison result (Fig. 2) showed that the normalized signals detected by these two systems were in good agreement, which suggested that B-mode imaging was as effective as the single-transducer ACD system. The experimental data were also simulated based on integrated calculation of the scattering cross-section. By setting three fitting parameters (viz., the mean radius of bubbles [R ₀], the standard deviation of bubble size distribution [SD], and the gas diffusion constant [D _A]), we assumed that the best simulation result would be obtained when the SDE achieved the minimum value. The results showed that the experimental results for both two systems could be best simulated with R ₀ = 15 µm, SD = 4 µm, and D _A = 2.5 × 10^–9 m²/s.

Although the results showed that the experimental results could be simulated very well with appropriate fitting parameters, it was of interest to determine if the simulation result was or was not unique, because there were so many parameters involved in the simulation. Figure 8 illustrates the dissolution curve detected by the B-mode system and two simulation curves obtained while using different parameter sets. The SDEs of these two simulation results were both within the range of (1 + 5%) × SDE_min. Obviously, similar fitting results could be achieved by introducing different parameter sets. For example, if a larger bubble (e.g., R ₀ = 26 µm) was assumed and the diffusion constant D _A increased correspondingly (e.g., 10 e^–9), the simulated decay curve looked approximately the same as the others. So, how can one determine which simulation result is correct? We should notice that, in Fig. 8, one of the two simulation curves assumed a much bigger bubble and much higher diffusion constant. The air diffusion constant in water under standard conditions (25  C, 1 atm, 1 atm = 1.01325 × 10⁵ Pa) is 2.4 × 10^–9 m²/s. Therefore, it is applicable to achieve best simulation by adopting the fitting parameters (e.g., R0 = 15 µm, SD = 4 µm, and D _A = 2.5 × 10^–9 m²/s) that can generate the minimum SDE. Thus, with the best fitting approach, we could at least get a relatively accurate estimation for the bubble size distribution and the diffusion constant. The parameter SD which presents the deviation of bubble size helps to get the better estimation, although it might not affect the trend of simulation.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. The dissolution curve measured by B-mode ultrasound. Similar simulation results can be obtained while using different parameter sets.

The above discussion suggests that it is possible to estimate the bubble size distribution and gas diffusion constant by simulating the experimental data properly. Thus, we would like to estimate the size distribution of the cavitation bubbles generated during in vivo SW studies by applying this method to the in vivo ESWL and ESWT data. Figure 3 shows that, for the ESWL experiment with 2-Hz PRF, the best fitting parameters were R ₀ = 35 µm, SD = 10 µm, and D _A = 5 × 10^–7 m²/s. At 0.5-Hz PRF, the fitting results obtained for the ESWL study (Fig. 4) were R ₀ = 20 µm, SD = 8 µm, and the D _A = 2.5 × 10^–7 µm²/s. Figure 7 illustrates that the best fitting parameters for the ESWT data were R ₀ = 12.5 µm, SD = 7 µm, and D _A = 3 × 10^–8 m²/s. Church^[28] and Sapozhnilov et al. ^[22] have modeled the dynamic behaviors of the SW-induced cavitation bubbles. Their calculation showed that SW pulses would drive these bubbles into a dramatic growth and collapse followed by a slow dissolution of the bubbles. They also reported that the calculated initial equilibrium radius of the dissolving cavitation bubble in water is ~ 30 µm, which was at the same order as our estimated bubble equilibrium radius. The maximum diffusion constant for O₂ saturated in whole blood could refer to the order of 10^–5. Because the dissolved O₂ concentration in kidney or capillaries could change with conditions, the estimated diffusion constants obtained here are considered to be in a reasonable range, which is greater than 10^–9 (air diffusion constant in water) and smaller than 10^–5 (O₂ diffusion constant in whole blood).

The best fitting results for the ESWL pig kidney studies (see Figs. 3 and 4) indicate that, at t = 0 s, larger equilibrium diameters of cavitation bubbles were detected by B-mode ultrasound at PRF = 2 Hz. This observation could result from bubble coalescence and is consistent with the results reported in previous literatures.^{[16,22,29,30]} First of all, we have to mention that, due to the slow frame rate of the B-model ultrasound, the R ₀ discussed in the current dissolution studies were the equilibrium sizes of residual bubbles after SW pulses. These residual bubbles might be formed by coalescence of small bubbles, which has been observed and described in the optical investigations of Postema et al. ^[29,30] They found that two bubbles will come into contact more quickly if their distance is short. Sapozhnikov et al. ^[22] reported that denser cavitation bubble clouds with larger radii could be observed with the use of higher PRFs. Tu et al. also reported similar phenomenon in the ESWL studies based on B-mode imaging quantification analyses.^[16] Therefore, it can be postulated that the bubble clouds will be denser when generated at higher PRF, which means more bubbles will be generated with shorter adjacent distance, which would promote bubble coalescence. As a result, the residual bubbles with larger equilibrium radius were generated eventually by the SWs with 2-Hz PRF.

Another interesting point we want to point out here is that, for the experimental data of in vivo ESWL at 2-Hz-PRF and ESWT, the tails of the experimental dissolution curves do not follow the simulation results that trend downward continuously. Previous researchers have reported that cavitation bubble nuclei can be stabilized in tissue against dissolution.^[31–34] The main possible mechanisms of bubble stabilization have been proposed as: (i) small spherical gas bubbles could be stabilized by a layer of surface-active materials; essentially shrinking bubbles might be stabilized by a continuous surface formed by molecules;^[33] and (ii) cavitation bubble nuclei could be prevented from dissolution if they form inside a crevice.^[34] The experimental data tends to approach some stationary minimum values finally, which implied that some forms of “bubble stabilization” might be involved in the dissolution procedure, especially close to the end of the process. However, the experimental data of in vivo ESWL at 0.5-Hz-PRF met the simulation result pretty well. As mentioned in the researches of Postema et al.,^[29,30] sound pressure could drive two bubbles into each other, furthermore, if the distance between two bubbles gets shorter, they may come into contact with each other faster. Meanwhile, larger and more densely populated cavitation bubble clouds were observed in SWL at higher PRF,^[22] which in turn had shorter distance between bubbles. By their studies, one could conclude that there were less cavitation events and little residual bubbles at slower PRF. Therefore, the dissolution process at slower PRF (e.g., 0.5 Hz) followed the simulation results better than higher PRF (e.g., 2 Hz).

This study investigated dissolution processes of cavitation bubbles generated during in vivo SW -induced treatments. Firstly, the feasibility of detecting the bubble dissolution process based on the B-mode imaging technique was confirmed by comparing with the ACD system in the in vitro experiments. Secondly, the SW-induced cavitation bubbles dissolution processes were detected by the B-mode imaging system during in vivo SW treatments (ESWL and ESWT). Finally, the dissolution processes were simulated to study the characteristics and behaviors of the residual bubbles. The results showed that (I) B-mode imaging technology is an effective tool to monitor the temporal evolution of cavitation bubbles dissolution procedures after the SW pulses ceased, which is important for evaluation and controlling the cavitation activity generated during subsequent SW treatments within a treatment period; (II) the cavitation bubble dissolution procedure detected by B-mode ultrasound can be simulated based on the calculation of the scattering cross-section with employing gas bubble dissolution equations and Gaussian bubble size distribution; (III) the characteristics of the bubbles, such as the bubble size distribution and gas diffusion, can be estimated by simulating the experimental data properly. In summary, this study showed the promising future of estimating the size distribution of SW-induced cavitation bubbles and other involved coefficients, by simulating the bubble dissolution procedures. This work should be helpful for understanding the time-varying SW-induced cavitation activity, which is important for controlling the effects of SW treatments. Furthermore, it may provide a possible method to achieve an optimal SW treatment of minimizing tissue injury without compromising the therapeutic effect by adjusting the dissolving process of residual cavitation bubbles between two sequent SWs.